Instalamos los paquetes necesarios
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(sf)
Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
library(gstat)
library(geoR)
'RandomFieldsUtils' will NOT use OMP
'RandomFields' will NOT use OMP
--------------------------------------------------------------
Analysis of Geostatistical Data
For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
geoR version 1.8-1 (built on 2020-02-08) is now loaded
--------------------------------------------------------------
library(spdep)
Loading required package: sp
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`
library(DataExplorer)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
library(psych)
library(ggcorrplot)
Loading required package: ggplot2
Attaching package: ‘ggplot2’
The following objects are masked from ‘package:psych’:
%+%, alpha
library(biogeo)
Attaching package: ‘biogeo’
The following object is masked from ‘package:spData’:
world
library(ggplot2)
library(leaflet)
Cargamos los datasets
estaciones <- read.csv("data_smn/preprocessed/estaciones_smn_v2.csv")
horarios <- read.csv("data_smn/preprocessed/datohorario20210420.csv")
Observamos como esta compuesto el dataset
glimpse(estaciones)
Rows: 123
Columns: 9
$ NOMBRE <chr> "BASE BELGRANO II", "BASE CARLINI (EX JUBANY)", "BASE ESPERANZA", "BASE MARAMBIO", "BASE ORC…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "BUENOS AIRES"…
$ LATITUD_GRADOS <int> 77, 62, 63, 64, 60, 68, 36, 38, 37, 36, 34, 38, 37, 36, 34, 34, 34, 34, 36, 37, 34, 34, 34, …
$ LATITUD_MINUTOS <int> 52, 14, 24, 14, 45, 8, 50, 44, 43, 12, 32, 1, 26, 21, 36, 49, 33, 58, 2, 56, 33, 41, 40, 27,…
$ LONGITUD_GRADOS <int> 34, 58, 56, 56, 44, 67, 59, 62, 59, 61, 58, 61, 61, 57, 58, 58, 60, 57, 59, 57, 58, 58, 58, …
$ LONGITUD_MINUTOS <int> 34, 38, 59, 43, 43, 8, 53, 10, 47, 4, 40, 20, 53, 44, 36, 32, 55, 54, 8, 35, 49, 44, 38, 53,…
$ ALTURA <int> 256, 11, 24, 198, 12, 7, 147, 83, 207, 94, 26, 247, 233, 9, 12, 20, 81, 23, 36, 21, 32, 36, …
$ NRO <int> 89034, 89053, 88963, 89055, 88968, 89066, 87641, 87750, 87649, 87640, 87570, 87683, 87637, 8…
$ NroOACI <chr> "SAYB", "SAYJ", "SAYE", "SAWB", "SAYO", "SAYS", "SAZA", "SAZB", "SAZJ", "SAZI", "SADO", "SAB…
glimpse(horarios)
Rows: 2,041
Columns: 8
$ FECHA <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20…
$ HORA <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5…
$ TEMP <dbl> 20.7, 19.8, 19.5, 19.2, 19.2, 18.9, 19.0, 18.8, 19.0, 19.1, 19.9, 20.3, 21.8, 22.9, 23.2, 23.7, 23.2, …
$ HUM <int> 70, 79, 79, 81, 81, 84, 84, 84, 84, 84, 81, 82, 68, 72, 68, 66, 71, 68, 75, 72, 73, 68, 68, 68, 90, 88…
$ PNM <dbl> 1020.8, 1020.9, 1020.9, 1020.4, 1020.0, 1020.3, 1020.7, 1021.1, 1021.6, 1022.2, 1022.4, 1021.9, 1021.7…
$ DD <int> 70, 70, 70, 70, 70, 50, 70, 50, 50, 50, 50, 50, 50, 70, 70, 70, 70, 70, 90, 90, 90, 90, 70, 70, 50, 50…
$ FF <chr> "28", "28", "26", "26", "22", "22", "17", "17", "15", "19", "20", "24", "20", "19", "20", "17", "15", …
$ NOMBRE <chr> "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPA…
Observamos que la variable que representa la velocidad del viento en km/h es de tipo Char, por lo cual debemos convertila en int.
horarios$FF <- as.numeric(horarios$FF)
NAs introduced by coercion
summary(estaciones$ALTURA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.0 53.5 145.0 331.5 455.0 3459.0
summary(horarios$HORA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 7.00 12.00 12.16 18.00 23.00
summary(horarios$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-20.00 16.10 20.20 19.09 23.80 33.60 1
summary(horarios$HUM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
18.00 59.00 72.00 71.18 85.00 100.00 4
summary(horarios$PNM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
960.8 1011.0 1015.1 1075.3 1019.2 3133.0 23
summary(horarios$DD)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 50.0 70.0 118.9 160.0 990.0 1
summary(horarios$FF)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 9.00 15.00 16.06 22.00 107.00 1
hist(estaciones$ALTURA, main = "Histograma de la altura", xlab = "Altura")
hist(horarios$HORA, main = "Histograma de horarios", xlab = "Horarios")
hist(horarios$TEMP, main = "Histograma de temperatura", xlab = "Temperatura")
hist(horarios$HUM, main = "Histograma de humedad", xlab = "Humedad")
hist(horarios$PNM, main = "Histograma de presion atmosferica", xlab = "Presion Atmosferica")
hist(horarios$DD, main = "Histograma de la direccion del viento", xlab = "Direccion del viento")
hist(horarios$FF, main = "Histograma de la velocidad del viento", xlab = "Velocidad del viento")
Creamos unos boxplot para visualizar la distribucionde datos que potencialmente nos interecen para proceder con el analisis
boxplot(estaciones$ALTURA, main = "Boxplot de la altura", xlab = "Altura")
boxplot(horarios$HORA, main = "Boxplot de horarios", xlab = "Horarios")
boxplot(horarios$TEMP, main = "Boxplot de temperatura", xlab = "Temperatura")
boxplot(horarios$HUM, main = "Boxplot de humedad", xlab = "Humedad")
boxplot(horarios$PNM, main = "Boxplot de presion atmosferica", xlab = "Presion Atmosferica")
boxplot(horarios$DD, main = "Boxplot de la direccion del viento", xlab = "Direccion del viento")
boxplot(horarios$FF, main = "Boxplot de la velocidad del viento", xlab = "Velocidad del viento")
#plot_boxplot(horarios, by = "HORA")
#plot_boxplot(horarios, by = "TEMP")
Ahora, veamos que tan correlacionadas estan estas variables
corr <- cor(horarios[, c(2,3,4,5,6,7)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Hay datos nulos en estos datasets?
print("Nulos en horarios")
[1] "Nulos en horarios"
which(is.na(horarios))
[1] 4853 6710 6894 7999 8000 8401 8935 9137 9138 9139 9140 9141 9142 9143 9144 9556 9557 9558 9559
[20] 9560 9561 9562 9563 9564 9565 9566 9567 9605 10976 13017
print("Nulos en estaciones")
[1] "Nulos en estaciones"
which(is.na(estaciones))
integer(0)
Borramos los nulos del dataset
estaciones = na.omit(estaciones)
horarios = na.omit(horarios)
Analizamos la simetria de la variable que representa la velocidad del vientoo
# Analisis de simetria
#library(psych)
skew(horarios$FF)
[1] 1.591713
kurtosi(horarios$FF)
[1] 7.23501
La medida de asimetria y kurtosi terminan de validar lo que observamos en el histograma. La variable FF es asimetrica a derecha y tiene una mayor concentracion de valores muy cerca de la media de la distribución y muy lejos de la cola de la distribucion.
Necesitamos convertir las latitudes y longitudes a un valor decimal para poder graficarlo en el mapa con leaflet
#library(biogeo)
latitud <- c()
longitud <- c()
for(i in 1:nrow(estaciones)) {
latitud[i] <- dms2dd(dd = estaciones[i, "LATITUD_GRADOS"], mm = estaciones[i, "LATITUD_MINUTOS"], ss = 0, ns = "S")
longitud[i] <- dms2dd(dd = estaciones[i, "LONGITUD_GRADOS"], mm = estaciones[i, "LONGITUD_MINUTOS"], ss = 0, ns = "S")
}
estaciones['LATITUD'] <- latitud
estaciones['LONGITUD'] <- longitud
Unimos las dos tablas usando la variable NOMBRE como punto para combinar los datasets
data <- full_join(estaciones, horarios, by = c("NOMBRE" = "NOMBRE"))
glimpse(data)
Rows: 2,040
Columns: 18
$ NOMBRE <chr> "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRA…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "…
$ LATITUD_GRADOS <int> 77, 77, 77, 77, 77, 77, 77, 77, 62, 62, 62, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 63, …
$ LATITUD_MINUTOS <int> 52, 52, 52, 52, 52, 52, 52, 52, 14, 14, 14, 14, 14, 14, 14, 14, 24, 24, 24, 24, 24, 24, 24, …
$ LONGITUD_GRADOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 58, 58, 58, 58, 58, 58, 58, 58, 56, 56, 56, 56, 56, 56, 56, …
$ LONGITUD_MINUTOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 38, 38, 38, 38, 38, 38, 38, 38, 59, 59, 59, 59, 59, 59, 59, …
$ ALTURA <int> 256, 256, 256, 256, 256, 256, 256, 256, 11, 11, 11, 11, 11, 11, 11, 11, 24, 24, 24, 24, 24, …
$ NRO <int> 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89053, 89053, 89053, 89053, 89053, 8…
$ NroOACI <chr> "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYJ", "SAYJ", "SAYJ", "SAY…
$ LATITUD <dbl> -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -62.…
$ LONGITUD <dbl> -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -58.…
$ FECHA <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20…
$ HORA <int> 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 1, 2,…
$ TEMP <dbl> -20.0, -16.2, -16.9, -13.0, -15.1, -16.6, -15.3, -16.8, -0.1, -0.5, -0.9, -0.5, 0.3, 1.1, 1.…
$ HUM <int> 70, 72, 80, 73, 86, 66, 77, 85, 88, 82, 77, 75, 84, 92, 99, 91, 43, 49, 67, 81, 67, 62, 77, …
$ PNM <dbl> 975.7, 970.4, 966.7, 964.2, 961.4, 963.0, 964.8, 967.1, 992.9, 993.0, 992.9, 994.6, 995.3, 9…
$ DD <int> 140, 140, 110, 140, 160, 160, 90, 110, 230, 230, 270, 290, 290, 320, 320, 320, 160, 200, 50,…
$ FF <dbl> 13, 30, 50, 52, 74, 61, 44, 13, 15, 15, 26, 28, 41, 37, 50, 74, 20, 7, 13, 7, 11, 46, 56, 74…
Eliminamos las observaciones nulas en el dataset
data = na.omit(data)
Volvemos a graficar el grafico de correlacion incluyendo la variable altura.
corr <- cor(data[, c(7,13,14,15,16,17,18)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Agrupamos los datos por nombre calculando el promedio y desvio del viento
data_agg = data %>%
group_by(NOMBRE) %>%
summarise(MEAN_VIENTO_KMH = mean(FF), SD_VIENTO_KMH = sd(FF), LONGITUD = unique(LONGITUD), LATITUD = unique(LATITUD), .groups = "keep")
Convertimos data_agg a data.frame ya que necesitamos este tipo de dato para poder trabajar
df_data_agg = data.frame(data_agg)
Transformamos df_data_agg en un archivo geográfico utilizando el código de proyección mercator
data_sf = st_as_sf(df_data_agg, coords = c("LONGITUD", "LATITUD"), crs = 4326)
as_Spatial(data_sf)
Error in as_Spatial(data_sf) : object 'data_sf' not found
Validamos la clase del nuevo dataframe
class(data_sf)
[1] "sf" "data.frame"
##TODO: graficar el mapa de densidades con los puntos de las estaciones en data.
En el siguiente grafico de la republica argentina se observan en color azul las estaciones meteorologicas donde se realizaron las mediciones de la variables
#library(leaflet)
# Queremos que en el mapa se vea como etiqueta el nombre de la base meteorologica. Para eso aplicamos la siguiente funcion:
labs <- lapply(seq(nrow(data_sf)), function(i) {
paste0( '<p>', data_sf[i, "NOMBRE"], '</p>' ) })
leaflet() %>%
addTiles() %>%
addCircles(data = data_sf, weight = 3, label = lapply(labs, htmltools::HTML))
NA
En el mapa observamos que hay puntos muy distantes de la argentina continental. Dado el proposito de este estudio, el cual es determinar la locacion geografica optima en base a la variable velocidad del viento, decidimos remover estas observaciones ya que no aporta informacion util y ademas agregan ruido a nuestro analisis.
Primero, vamos a borrar las estaciones que no estan en la argentina continental - Base Carlini - Base San Martin - Base Marambio - Base Esperanza - Base Orcadas
data_sf = data_sf[data_sf$NOMBRE != "BASE CARLINI (EX JUBANY)",]
data_sf = data_sf[data_sf$NOMBRE != "BASE SAN MARTIN",]
data_sf = data_sf[data_sf$NOMBRE != "BASE MARAMBIO",]
data_sf = data_sf[data_sf$NOMBRE != "BASE ESPERANZA",]
data_sf = data_sf[data_sf$NOMBRE != "BASE ORCADAS",]
data_sf = data_sf[data_sf$NOMBRE != "BASE BELGRANO II",]
Tambien vimos que El Bolson tiene las 9 observaciones en 0 para la variable FF y DD. Consideramos esto como un error en el instrumento de medicion, por lo cual vamos a eliminar a esa estacion del analisis.
data_sf = data_sf[data_sf$NOMBRE != "EL BOLSON AERO",]
Repetimos el plot para validar
labs <- lapply(seq(nrow(data_sf)), function(i) {
paste0( '<p>', data_sf[i, "NOMBRE"], '</p>' ) })
leaflet() %>%
addTiles() %>%
addCircles(data = data_sf, weight = 3, label = lapply(labs, htmltools::HTML))
NA
Analisis univariado
describe(data_sf$MEAN_VIENTO_KMH)
hist(data_sf$MEAN_VIENTO_KMH)
boxplot(data_sf$MEAN_VIENTO_KMH)
hist(data_sf$SD_VIENTO_KMH)
boxplot(data_sf$SD_VIENTO_KMH)
#COMMENTS: sacar inliers (manual de buenas practicas compartido por silvia)
hist(data_sf$MEAN_VIENTO_KMH)
boxplot(data_sf$MEAN_VIENTO_KMH)
qqnorm(data_sf$MEAN_VIENTO_KMH)
qqline(data_sf$MEAN_VIENTO_KMH, col=2)
#TODO: sacar outliers y puntos en la arg continental.
Hacemos el test de normalidad contra la variable MEAN_VIENTO_KHM
shapiro.test(data_sf$MEAN_VIENTO_KMH)
Shapiro-Wilk normality test
data: data_sf$MEAN_VIENTO_KMH
W = 0.9698, p-value = 0.03265
El p-value obtenido es menor a 0.05, lo cual indica que los datos que tenemos no son normales
knea <- knearneigh(data_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran
Moran I test under randomisation
data: data_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 4.0432, p-value = 2.636e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.52033029 -0.01111111 0.01727672
moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)
data_sf[10,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -65.76667 ymin: -28.6 xmax: -65.76667 ymax: -28.6
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
16 CATAMARCA AERO 46.28571 11.51397 POINT (-65.76667 -28.6)
data_sf[56,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -62.38333 ymin: -37.6 xmax: -62.38333 ymax: -37.6
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
63 PIGUE AERO 38.23077 6.905962 POINT (-62.38333 -37.6)
data_sf[87,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -61.96667 ymin: -33.66667 xmax: -61.96667 ymax: -33.66667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
94 VENADO TUERTO AERO 32 11.09955 POINT (-61.96667 -33.66667)
data_sf[21,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -61.88333 ymin: -37.43333 xmax: -61.88333 ymax: -37.43333
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
27 CORONEL SUAREZ AERO 27.33333 5.278889 POINT (-61.88333 -37.43333)
data_sf[73,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -66.35 ymin: -33.26667 xmax: -66.35 ymax: -33.26667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
80 SAN LUIS AERO 17 7.901568 POINT (-66.35 -33.26667)
data_sf[51,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -64.31667 ymin: -23.15 xmax: -64.31667 ymax: -23.15
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
58 ORAN AERO 1.222222 2.981424 POINT (-64.31667 -23.15)
data_sf[55,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -71.01667 ymin: -46.51667 xmax: -71.01667 ymax: -46.51667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
62 PERITO MORENO AERO 1.1 3.478505 POINT (-71.01667 -46.51667)
data_sf = data_sf[data_sf$NOMBRE != "CERES AERO",]
data_sf = data_sf[data_sf$NOMBRE != "SAN LUIS AERO",]
data_sf = data_sf[data_sf$NOMBRE != "PERITO MORENO AERO",]
data_sf = data_sf[data_sf$NOMBRE != "ORAN AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CORONEL SUAREZ AERO",]
data_sf = data_sf[data_sf$NOMBRE != "PIGUE AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CATAMARCA AERO",]
# Creamos una lista de vecinos
knea <- knearneigh(data_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
# Hacemos el test de moran
globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran
Moran I test under randomisation
data: data_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 4.8984, p-value = 4.83e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.66907198 -0.01204819 0.01933449
moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)
data_sf[81,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -63.01667 ymin: -40.85 xmax: -63.01667 ymax: -40.85
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
95 VIEDMA AERO 16.5 5.890302 POINT (-63.01667 -40.85)
data_sf[60,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -69.28333 ymin: -51.61667 xmax: -69.28333 ymax: -51.61667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
73 RIO GALLEGOS AERO 25.91667 6.907944 POINT (-69.28333 -51.61667)
data_sf[80,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -61.96667 ymin: -33.66667 xmax: -61.96667 ymax: -33.66667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
94 VENADO TUERTO AERO 32 11.09955 POINT (-61.96667 -33.66667)
data_sf[19,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -58.76667 ymin: -27.45 xmax: -58.76667 ymax: -27.45
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
28 CORRIENTES AERO 9.916667 5.115336 POINT (-58.76667 -27.45)
shapiro.test(data_sf$MEAN_VIENTO_KMH)
Shapiro-Wilk normality test
data: data_sf$MEAN_VIENTO_KMH
W = 0.96946, p-value = 0.04325
data_sf = data_sf[data_sf$NOMBRE != "VIEDMA AERO",]
data_sf = data_sf[data_sf$NOMBRE != "RIO GALLEGOS AERO",]
data_sf = data_sf[data_sf$NOMBRE != "VENADO TUERTO AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CORRIENTES AERO",]
# Creamos una lista de vecinos
knea <- knearneigh(data_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
# Hacemos el test de moran
globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran
Moran I test under randomisation
data: data_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 4.8722, p-value = 5.518e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.67737208 -0.01265823 0.02005783
moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)
plot(data_sf$MEAN_VIENTO_KMH)
plot(log(data_sf$MEAN_VIENTO_KMH))
qqnorm(log(data_sf$MEAN_VIENTO_KMH))
qqline(log(data_sf$MEAN_VIENTO_KMH), col=2)
shapiro.test(data_sf$MEAN_VIENTO_KMH)
Shapiro-Wilk normality test
data: data_sf$MEAN_VIENTO_KMH
W = 0.97049, p-value = 0.06096
#coordinates(data_sf) = ~LATITUD+LONGITUD
#class(data_sf)
plot(df_data_agg)
plot(data_agg[, c(4,5)])
v <- variogram(MEAN_VIENTO_KMH~1, data_sf)
plot(v)
v <- variogram(MEAN_VIENTO_KMH~1, data_sf)
plot(v)
vt_exp = fit.variogram(v, vgm(125, "Exp", 30, 5))
vt_exp
plot(v , vt_exp)
vt_mat = fit.variogram(v, vgm(125, "Mat", 30, 5))
plot(v , vt_mat)
vt_exc = fit.variogram(v, vgm(125, "Exc", 30, 5))
No convergence after 200 iterations: try different initial values?
plot(v , vt_exc)
Kriging
attr(vt_exp, 'SSErr')
[1] 0.1189922
attr(vt_mat, 'SSErr')
[1] 0.1189922
attr(vt_exc, 'SSErr')
[1] 11.92923
attr(vt_bes, 'SSErr')
[1] 0.1344978
ko1 <- krige(MEAN_VIENTO_KMH~1, data_sf, grilla, model = vt_exp, nmax=20)
spplot(ko1["var1.pred"])
spplot(ko1["var1.var"])